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A  VISCOELASTIC  BEM  FOR  MODELING  OXIDATION 


Thye-Lai  Tung,  Jerome  Connor,  and  Dimitri  A.  Antoniadis 
Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts  02139 

Abstract 

A  viscoelastic  boundary  element  method  has  been  developed  to  model  the  motion 
of  silicon  dioxide  and  silicon  nitride  during  thermal  oxidation  of  silicon.  This  tech¬ 
nique  uses  Kelvin’s  solution  reformulated  according  to  the  correspondence  principle 
on  viscoelasticity.  Constant-velocity  loading  is  chosen  to  ensure  smooth  variations  in 
displacement  and  stress  behavior  for  a  wide  range  of  relaxation  times. 

1  Introduction 

A  major  problem  in  modeling  thermal  oxidation  of  silicon  is  the  moving  boundaries.  Driven 
by  the  conversion  of  silicon  to  silicon  dioxide  of  increased  volume,  boundaries  change  shape 
drastically  during  a  local  oxidation  of  silicon  (LOCOS)  process  step.  Conventional  methods 
such  as  the  finite  element  method  (FEM)  require  a  mesh  to  subdivide  the  simulation 
domain.  To  apply  such  methods  to  thermal  oxidation,  one  must  first  develop  computer 
codes  to  regenerate  the  mesh  automatically  and  optimally  as  the  oxide  changes  shape  at 
every  time  step  [1,2,3].  After  that,  one  may  have  to  deal  with  the  transfer  of  stress  history 
from  the  old  mesh  to  the  new  one,  depending  on  what  oxide  flow  model  is  used.  This  is 
another  difficult  problem  if  one  wishes  to  minimize  numerical  diffusion  of  stress  distribution 
due  to  regriding. 

In  contrast  to  the  FEM,  the  boundary  element  method  (BEM)  does  not  require  a 
mesh  in  general  since  it  does  all  calculations  on  the  boundary  but  not  in  the  interior.  Be¬ 
cause  segments  need  not  be  regenerated,  keeping  track  of  stress  history  on  the  boundary 
is  straight  forward.  However,  previous  efforts  on  applying  the  BEM  to  model  oxide  mo¬ 
tion  as  viscoelastic  flow  suffered  from  some  drawbacks.  In  Matsumoto’s  formulation  [4], 
domain  calculations  were  needed.  Although  restrictions  on  the  grid  were  not  as  stringent 
as  those  in  the  FEM,  excessive  computation  time  was  evident.  Simulation  setup  might 
still  be  difficult  for  complex  structures.  Isomae  attacked  the  problem  differently  by  using 
a  Laplace  transform  technique  [5|.  According  to  the  correspondence  principle  of  linear 
viscoelasticity,  we  may  solve  a  viscoelastic  problem  through  an  equivalent  elastostatic  sys¬ 
tem  in  the  Laplace  transform  space.  Boundary  element  techniques  for  elastostatics  do  not 
require  domain  calculations,  but  unfortunately  it  is  impossible  to  solve  Laplace  transforms 
analytically  in  any  practical  simulation  situations.  Numerical  Laplace  transform,  as  used 
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by  Isomae,  is  deemed  wasteful  because  a  new  solution  must  be  generated  at  every  time 
step. 

In  this  paper  we  describe  a  generalized  BEM  for  two-dimensional  linear  viscoelastic 
flow,  with  application  to  thermal  oxidation.  Instead  of  performing  the  Laplace  transform 
on  the  global  solution,  we  operate  it  on  the  fundamental  equations.  Our  kernel  functions 
are  derived  from  Kelvin’s  solution  used  in  elastostatic  BEM  applications  [6].  A  suitable 
excitation  function  is  added  to  produce  a  time-varying  body  force  that  causes  the  load 
point  to  move  with  a  constant  velocity  in  a  viscoelastic  medium.  This  approach  eliminates 
domain  calculations  and  numerical  Laplace  transforms.  Yet,  it  is  as  easy  to  use  as  those 
for  elastostatics.  The  main  disadvantage  is  that  boundary  conditions  cannot  be  satisfied 
for  all  times,  so  the  system  must  be  solved  periodically.  But  that  is  not  a  problem  for 
thermal  oxidation  because  a  new  solution  must  be  obtained  anyway  for  every  time  step. 
This  formulation  can  handle  all  possible  values  of  Poisson’s  ratio  and  a  wide  range  of 
stress  relaxation  times,  essentially  encompassing  viscous  incompressible  flow  and  elastic 
deformation. 

2  Numerical  Formulation 

We  implement  the  viscoelastic  BEM  using  the  so-called  indirect  formulation,  which  is 
also  known  as  the  classical  or  source  method.  Here  we  attempt  to  model  a  problem 
by  putting  “sources”  on  the  boundary  and  adjust  their  strength  so  that  the  fields  they 
generate  match  the  prescribed  boundary  conditions.  For  a  potential  problem,  the  sources 
are  simply  electrical  charges.  The  indirect  formulation  has  an  advantage  over  the  direct 
formulation  in  that  different  components  within  the  stress  tensor  or  displacement  vector 
can  be  obtained  more  readily. 

The  thermal  oxidation  process  actually  consists  of  two  tightly  coupled  processes,  namely 
oxidant  diffusion,  and  oxide  motion.  For  details  on  the  oxidant  diffusion  process,  boundary 
conditions,  and  numerical  implementation  of  the  integral  equations,  readers  are  referred 
to  our  previous  paper  on  the  same  subject  [7].  In  that  paper  we  use  an  incompressible 
viscous  flow  model  for  oxide  motion;  here  we  have  a  generalized  model  applicable  to  the 
oxide,  silicon  nitride  and  silicon  substrate. 

The  two  sets  of  integral  equations  for  viscoelasticity  are  shown  below: 


■“(*)  =  e)dr 


where  (r,-y,  uy,  and  p  are  the  stress  tensor,  displacement  vector,  and  source  density  respec¬ 
tively.  r  denotes  the  boundary.  and  u,*y  are  the  kernels;  their  actual  forms  are  given 
in  the  appendix.  For  compactness,  Einstein  notation  has  been  used.  To  get  the  surface 
traction,  we  apply  the  following  formula: 

Pi{^]  = 

where  ny  is  the  direction  cosine  of  the  surface. 

If  we  want  to  treat  silicon  nitride  as  an  elastic  material  with  stiffness,  we  have  to  modify 
the  boundary  condition  of  the  oxide-nitride  interface  from  a  free  surface  condition  p  =  0 
to  the  following: 

1.  displacement  vector  is  continuous  across  the  interface, 

2.  surface  tractions  of  the  two  materials  are  equal  and  opposite. 

To  ensure  stability,  the  equations  for  silicon  nitride  layer  and  oxide  region  must  be  solved 
simultaneously.  The  resulting  system  matrix  is  larger  but  banded. 

Likewise,  we  use  the  same  approach  if  we  desire  to  model  the  silicon  substrate  as  an 
elastic  foundation,  instead  of  treating  it  as  a  rigid  body.  Note  that  due  to  the  unique 
formulation  of  BEM,  boundary  conditions  at  y  =  — oo  need  not  be  specified  for  the  silicon 
substrate. 

In  viscoelasticity,  we  must  include  past  stress  history  in  the  present  time  step.  Consider 
p”,  the  normal  component  of  the  surface  traction  at  time  step  m.  The  two  different  relax¬ 
ation  terms,  as  defined  in  the  appendix,  are  kept  seperated  and  updated  in  the  following 
way: 

Pna  =  Pna  +  Pna"*  eXp(-Ar/r„) 

Pn^  =  Pn0+Pni^eXp{-AriTp) 

where  At”*  is  the  time  step  size,  Ta  and  Tp  the  characteristic  relaxation  time  constants. 

3  Simulation  Results 

We  will  demonstrate  how  stress  is  relieved  via  viscoelastic  flow,  the  effect  of  the  silicon 
nitride  layer  on  the  oxide  shape  and  the  stress  distribution.  The  relaxation  time,  or 
equivalently  the  viscosity,  of  oxide  has  not  been  determined  accurately.  It  is  diffucult  to 
'  do  so  because,  for  a  given  temperature,  the  viscosity  depends  on  the  quality  of  the  oxide, 
water  contents,  and  the  presence  of  other  impurities.  Most  of  the  recent  data  on  oxide 


viscosity  are  inferred  rather  than  measured  directly  in  experiments,  but  they  are  useful  as 
ballpark  figures.  In  any  case,  it  is  a  good  exercise  to  vary  the  relaxation  time  to  see  how 
it  affects  the  stress  distribution. 

Shown  in  Fig.  1  is  the  outline  of  a  semi-recessed  LOCOS  structure  plotted  for  every 
time  step.  In  this  simulation  window  of  1.6/im  wide,  the  silicon  nitride  mask  extends  from 
X  =  0  to  X  =  0.96/xm  and  is  assumed  to  be  totally  flexible.  Initially  the  structure  has  a  pad 
oxide  thickness  of  200A.  Oxidation  is  carried  out  at  925‘‘C  in  a  wet  ambient  for  3.4  hours 
to  get  a  final  field  oxide  thickness  of  SOOOA.  The  Young’s  modulus  and  Poisson’s  ratio  are 
taken  to  be  8  x  10^^  dynes-cm~’  and  0.194.  Shown  in  Fig.  2a,  2b  and  2c  is  the  normal 
surface  traction  at  the  oxide-silicon  interface  corresponding  to  relaxation  times  (^)  of  100 
hours,  1  hour,  and  1  minute  respectively.  The  normal  component  of  the  surface  traction  is 
plotted  for  every  time  step,  just  like  the  outline  of  the  oxide  in  Fig.  1.  (Because  the  oxide 
shape  is  almost  identical  for  all  relaxation  times,  only  one  is  shown  in  Fig.  1.)  In  all  the  3 
stress  plots,  there  are  two  peaks  in  the  compressible  stress  region  (negative  value  range). 
The  early  peak  occurs  at  the  edge  of  the  nitride  mask  (x  =  0.96/xm).  This  peak  is  due  to 
the  highly  nonuniform  oxidation  rate  in  that  region.  As  time  progresses,  the  peak  shifts 
to  the  left,  further  into  the  nitride  mask,  and  gives  rise  to  a  late  peak.  As  expected,  stress 
decreases  as  the  ’'iscosity  gets  lower. 

In  the  second  part,  we  repeat  the  same  simulations  with  the  silicon  nitride  layer  modeled 
as  an  elastic  material.  The  Young’s  modulus  and  the  Poisson’s  ratio  of  silicon  nitride  are 
assumed  to  be  3.29  x  10“  dynes-cm~*  and  0.266  respectively  |6).  The  thickness  of  the 
nitride  layer  is  0.1/xm.  The  final  shapes  of  the  oxide  and  the  nitride  layer  are  shown  in 
Fig.  3a,  3b,  and  3c.  As  we  can  see,  the  nitride  layer  bends  less  as  the  oxide  becomes  less 
viscous  and  flows  more  readily.  The  corresponding  values  for  peak  stress  are  2  x  10“, 
1.4  X  10“,  and  4.6  x  10‘°  dynes  cm”*  respectively.  The  first  two  values  are  unrealistically 
high  because  they  are  the  same  order  of  magnitude  as  the  elastic  moduli  of  oxide.  We  find 
that  in  general  we  need  those  stress  values  in  order  to  bend  the  nitride  mask  to  a  degree 
comparable  to  what  is  found  in  experiments.  To  keep  stress  down  to  a  realistic  level,  silicon 
nitride  must  deform  elastoplastically  or  viscoelastically.  Unfortunately  we  don’t  have  data 
on  those  behaviors. 

4  Summary 

A  boundary  element  technique  has  been  developed  to  model  thermal  oxidation  of  silicon 
in  two  dimensions.  It  can  handle  a  wide  range  of  relaxtion  times  for  the  viscoelsistic  flow 
of  oxide.  Simulations  of  some  simple  structures  have  been  demonstrated.  Preliminary 
results  indicate  that  it  is  inadequate  to  treat  silicon  nitride  as  an  elastic  material.  A 
comprehensive  characterization  of  silicon  nitride  is  clearly  needed. 
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Appendix 

The  fundamental  solutions  for  viscoelastic  flow  induced  by  constant-velocity  loading  are 
given  below.  Note  that  the  subscripts  of  ^  denote  partial  derivatives. 
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where  E  is  the  Young’s  modulus,  G  the  shear  modulus,  and  u  the  Poisson’s  ratio.  The 
fundamental  solutions  are  functions  of  time;  they  are  evaluated  at  the  end  of  a  time  step 
of  size  At.  Ka  and  Kp  decay  exponentially  with  a  time  constant  of  Ta  and  Tp  respectively 
when  the  load  is  removed. 
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Fig.  1.  Outline  of  oxide  for  every  time  step.  Final  field  oxide  thickness  is  0.5/xm. 
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